On the Effect of Constraint Enforcement on the Quality of 
Numerical Solutions in General Relativity 
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In Brodbeck et al 1999 it has been shown that the linearised time evolution 
equations of general relativity can be extended to a system whose solutions asymp- 
totically approach solutions of the constraints. 

In this paper we extend the non-linear equations in similar ways and investigate 
the effect of various possibilities by numerical means. Although we were not able 
to make the constraint submanifold an attractor for all solutions of the extended 
system, we were able to significantly reduce the growth of the numerical violation 
of the constraints. Contrary to our expectations this improvement did not imply a 
numerical solution closer to the exact solution, and therefore did not improve the 
(3JT)! quality of the numerical solution. 



I. INTRODUCTION 

Many physical theories are based upon systems of partial differential equations which 
contain more equations than variables like Maxwell's equations or general relativity. The 
initial data for the time evolution equations cannot be given freely, they must satisfy con- 
straints. It is necessary for the consistency of the theory, that for any data of the time 
evolution equations which initially satisfy the constraints, the constraints are satisfied for 
all times. This property is called "propagation of constraints" . 

Let us consider Maxwell's equations in vacuum as a simple example. The time evolution 
equations tell us that the time derivative of the electric and magnetic field are proportional 
to the curl of the magnetic and electric field. The vanishing of the divergence of the electric 
and magnetic field are the constraints. It can easily be shown that the constraints propagate. 
In cases where the solutions of a system of partial differential equations are determined by 
numerical means we cannot expect to get an exact propagation of the constraints. Due to the 
discretisation of the equations the numerical solution deviates from the exact solution by the 
discretisation error. As a consequence, the constraints are not fulfilled exactly after having 
evolved for some time, even if the initial data solved the constraints. In the spirit of M we 
call a discretisation of the time evolution equations compatible with the constraints, if the 
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numerical violation of the constraints has the same convergence order as the discretisation 
of the time evolution equations. Unfortunately the experience of numerical relativity shows 
that compatibility is not sufficient for obtaining numerical solutions with small numerical 
violations of the constraints. In many cases the violation of the constraints seems to grow 
at least exponentially with time. This effect is believed to be a major contribution to the 
numerical error of numerically calculated solutions. 

In this work we examine the effect of changing the evolution equations outside the sub- 
manifold of data on which the constraints are satisfied. Although our method is different, we 
should mention that already in || such a change has been suggested. Furthermore, to our 
knowledge, we perform the first systematic analysis of the correlation between the violation 
of the constraints and the quality of numerical solutions in general relativity. 
As the solutions of the field equations of general relativity satisfy the constraint equations for 
all times, the solutions are not affected by modifications of the evolution equations for data 
which do not satisfy the constraints. Let us denote the subspace of the function space of so- 
lutions to the evolution equations which satisfy the constraints as "constraint submanifold" . 
In |]J h has been proven, that at least for the linearised Einstein equations the constraint 
submanifold can be made an attractor for the linearised time evolution equations. 
If the solution of the evolution equations automatically approaches the constraint subman- 
ifold, the system of evolution equations carries a dissipative term in it, and therefore, the 
numerical solution will also approach the constraint submanifold provided the grid is not 
too coarse. Therefore, to avoid a numerical violation of the constraints, it is sufficient to 
make the constraint submanifold an attractor of the (modified) evolution equations. In such 
a case the constraint submanifold is 'asymptotically stable'. 

In Brodbeck et al. [[TJ a general method has been proposed to derive symmetric hyperbolic 
extensions of symmetric hyperbolic evolution equations with first order constraints which are 
promising candidates for asymptotic stability. These extended systems are called A-systems. 
In the same article it has also been proven, that at least in the case of the linearised Einstein 
equations there exist parameters such that the constraint submanifold is indeed an attractor 
for the modified evolution equations. 

As the extension of the analysis to the non-linear Einstein equations seemed to be beyond the 
scope of present analytical techniques, we took a numerical approach in this paper and in- 
vestigated the following questions: Firstly, can we, simply by way of numerical experiments, 
find a A-system for the non-linear Einstein equations for which the constraint submanifold 
is attractive? And secondly, is the numerical solution of the modified systems closer to 
the exact solution than the numerical solution of the unmodified system? To reduce the 
complexity and to have exact solutions available to compare with, we have restricted our 
investigations to solutions with two Killing vectors. 

In our experiments we were able to find a variety of A-systems for which the violation of all 
constraints is improved. However, we did not find a single system for which the constraint 
submanifold is asymptotically stable. Surprisingly, the improvement in the constraint vio- 
lation did not imply an improvement of the numerical solution. 

It is important to mention that a general attractive force towards the constraint submanifold 
does not guarantee the numerical solution to approach the exact solution corresponding to 
the initial data used. Regardless of the system of the field equations of general relativity, 
there are additional degrees of freedom which can be affected by the additional terms in the 
A-system. In our numerical experiments, these additional degrees of freedom were affected 
in such a way that in general the numerical solution of the modified system was not closer 
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to the exact solution, even if it was closer to the constraint submanifold. 

This paper is structured as follows: In chapter II we introduce parametrised A-systems 
and describe the simplifications implied by the symmetry assumptions. In the next chapter 
we sketch the numerical implementation, recall important features of the exact solutions 
used in the comparisons, and define the measures used for quality assessments. Chapter [IV] 
contains the actual numerical investigations, where we describe the performed probing of the 
parameter space. Using selected examples we demonstrate the influence of the individual 
parameters on the quality of the numerical solution. 



II. THE PARAMETRISATION OF THE A-SYSTEM 

The construction of a A-system is based on a split of the system of equations into sym- 
metric hyperbolic evolution equations and first order constraints []TJ. There are various 
possibilities to write Einstein's equations in a form like this. 

In our work, we take the conformal field equations ||, fj|, f| in the first-order formulation 
described in 0. Looking at the equations (13) and (14) of 0, it is easy to see that Einstein's 
equations and their extension, the conformal field equations, are a "quasilinear version" of 
Maxwell's equations. We use the conformal field equations instead of Einstein's equations 
to obtain an easy and well defined treatment of grid boundaries, as discussed in 0. Since 
we are primarily interested in the effect of the non-linearities, we can reduce the compu- 
tational complexity by restricting ourselves to asymptotically A3 spacetimes ||, which are 
spacetimes with two commuting, hypersurface orthogonal Killing vector fields. We align 
the y and z coordinates with the Killing orbits. Our solutions do therefore not depend on 
the spacelike coordinates y and z. Under these symmetry assumptions, the conformal field 
equations can be written in the following form 

A l 9 + A l 9 - 6 » = (la) 

|/-»/ = (lb) 

£/-«* = <>, do) 

with a time coordinate t labelling the spacelike hypersurfaces, x being the non-Killing space- 
like coordinate, and 

g = (kn, 7 1 n, E 22 , E 33 , B 23 , ™R 1} ™R U ) (2a) 
/ = (/in, h 22 , h 33 , k 22 , k 33 , 7S2, 7*33, E u , (M) i? 22 , (1 ' 1) i?33, 0, ^o, ^1,00). (2b) 

The tensor h a b is the 3-metric, k a b the extrinsic curvature of the spacelike hypersurfaces, 
7 a hc the connection for h a b, (0,1) ^a and (1,1) R a b parts of the tracefree part of the Ricci tensor, 
E a b and B a b the electric and magnetic part of the conformal Weyl tensor, Q the conformal 
factor, Q and Qi its normal and space derivative, and u a second derivative of the conformal 
factor, as described in more detail in [[?]]. The symmetric matrices A, A and the vectors b g , 
bj and Cf depend on g, f and gauge functions. The matrix A is positive definite, hence the 
system consisting of the equations (|Ta|) and (|Tb|) is symmetric hyperbolic. 
The variables, which are functions of t and x only, have been split into two classes, called 
g and /. For the variables denoted by / the system contains evolution equations and 
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constraints. Since there are only evolution equations for the variables g, these represent the 
degrees of freedom. 

When we evolve initial data forward in time by means of the evolution equations flla|) 
and flTED, the constraints ( |I~c| ) fulfil an evolution equation from which the propagation of the 



constraints can be derived. 

In the following we will call ( |Ta]) and flTED the 'unextended system'. To obtain the new, 
extended system, the A-system, additional variables, called A, are introduced and the con- 
straint equations are extended to evolution equations for the new variables: 

A^S+ b 3 + C\ = (3a) 

| /+ B±X- 6/ + DA = o (3b) 

I a+bT J^- bTc ' +ea = o - < 3c > 

The quantities B, C, D, and E are matrices. B T denotes the transposed matrix of B. This 
system is constructed in such a way that 

1. it is symmetric hyperbolic. 

2. in the case in which the variables A vanish identically the system is reduced to the 
original system ([[]). 

Due to the second requirement, the A-system is a generalisation of the original system. The 
first requirement implies well-posedness of the initial value problem. Apart from the restric- 
tions resulting from the two conditions above, the choice of the parameters B, C, D, and E 
is free. It is the aim to choose them in such a way that for all solutions of the system the 
variables A decay, which then implies that the solution is driven towards a solution of the 
constraints. 

We will now shortly explain why we have introduced these parameters. Let us assume that 
the constraint equations are not fulfilled exactly, i.e. d x f — Cf ^ 0, and A = initially. In the 
case of vanishing E, the variables A are the time integral of the violation of constraints - as 
a result of the new evolution equation Q5c|). For non- vanishing E, in addition, the EA-term 
will have a damping or amplifying effect onto the variables A, depending on the eigenvalues 
of E. For positive eigenvalues we expect damping, for negative eigenvalues amplification. 
The information about the violation of constraints, saved in the variables A, is coupled back 
to the variables (g, f) by the terms Hd x \, DA, and CA. 

As the conformal field equations are a generalisation of Einstein's equations we can relate 
the constraints flic]) , or more general, the constraints of the conformal field equations with- 
out symmetries 0, to the momentum and Hamiltonian constraints of the standard 3+1 
equations. The Hamiltonian constraint and the momentum constraint form a subset of the 
constraints used in our system. Contracting equation (14b) of [0] with the 3- metric h bc and 
restricting it to the case of Q — 1, we recover the momentum constraint (3) V&(/c ab — h ab k) = 0, 
with k = k a bh ab . Similarly, contracting equation (14c) of twice and evaluating the 3- 
Christoffel symbols, one can deduce the Hamiltonian constraint {3) R + k 2 — k a bk ab = 0, where 
{3) R denotes the Ricci scalar of the 3-metric. 
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III. THE NUMERICAL IMPLEMENTATION 

In this section we will describe the basic elements of the numerical implementation which 
we use to compare the quality of solutions obtained by solving the A-system with the quality 
of solutions obtained by solving the unextended system. These elements are the construction 
of initial data, the scheme to numerically integrate the time evolution equations, and the 
measures used to assess the quality of our numerical solutions. In addition we will briefly 
describe the exact solutions which we used as reference solutions. 



A. Constructing hyperboloidal initial data 

In order to analyse the numerical behaviour of the A-system, we first have to construct 
initial data for the conformal field equations. These data are called hyperboloidal initial 
data. 

In [||, section 2] the procedure of calculating initial data has been described in detail for the 
case without any symmetry assumptions. We slightly modified the procedure by making use 
of the symmetry assumptions, namely that our spatial grid is only one- dimensional (ID). 
For an exact solution we prescribe the 4-metric g a f, and the conformal factor Q as functions 
of (t,x). From those we calculate our variables (g,f) and the gauge source functions nu- 
merically. The code has also got the functionality to perform a coordinate transformation 
to express the exact solution in new coordinates (t',x f ). 

In the calculations presented in this paper we used A = as initial setting for A. 



B. The integration of the time evolution equations 

In order to discretise the evolution equation 

d t u + £(u) d x ii = b(u) (4) 

for the vector of variables u, we adjust the second-order scheme described in |5|, section 3.1] 
to our symmetry assumption, i.e. we perform a Strang splitting ansatz to split equation (ffl) 
into a principal part 

d t u + A{u) d x u = (5) 

and a source part 

d t u = b(u). (6) 
We then solve the principal part (|5|) by the rotated Richtmyer scheme. With one spatial 



dimension, this is equivalent to the standard second-order Lax-Wendroff method [10]. The 
source part is integrated by the pseudo-implicit Heun-scheme. As described in ||, principal 
and source part are combined in different order, depending on whether the time step is odd 
or even, to achieve second-order convergence. 

In || it was shown how superior a 4th order scheme would be. In a normal application 
this superiority would be a big advantage — here it is a disadvantage, however: We are going 
to analyse the impact of a drift away from the constraint submanifold on the quality of the 
solution. This drift originates in the discretisation error. If the scheme is very accurate, the 
drift is very small, too. But then, the differences in the quality of the numerical solutions 
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are also very small, which makes it harder to distinguish between them. This explains why 
we use the second-order scheme. 

In the runs described in part [TV] we cover the spacelike hypersurface with 161 grid points. 
The length of the time step is chosen dynamically by evaluating the Courant-Friedrichs- 
Lewy condition for each time slice. If not explicitly stated otherwise, we use half of what 
would be allowed by the Courant-Friedrichs-Lewy condition. 



C. Measures of quality 

We use the following measures to analyse the quality of a numerical solution. 
Firstly, we determine the numerical violation of the constraints as a measure of the distance 
from the constraint submanifold. To be able to present our findings with a limited number 
of plots we condense the information by using the "norm" 

l|A C ||(*) :=</£/ Ct(t,x)*dx, (7) 

where the summation includes all constraints Ci(t,x). In the actual numerical calculations 
the integral is of course replaced by a sum over all grid points with Q > 0, which represent 
physical spacetime. In the conformal approach to numerical relativity grid points with Q < 
represent a formal extension of the grid without physical relevance. Therefore, it would be 
wrong to include them into the measure. 

Secondly, we compare our numerical solution to the exact solution. Since the numerical cal- 
culation of ( (11> /2 a fe, E ao , B a b) from the given solution (<? a &, fl) involves solving elliptic equa- 
tions (cf. 0), we compare (fl (1,1) i? a &, Q 2 E a b, QB a t,) to the corresponding quantities from the 
exact solution and call the result "pseudodifference" Vi(t,x). In || it was found that with 
respect to the relative error this is equivalent to the difference in the variables. Again, to 
be able to present our findings with a limited number of plots, we condense the information 
by using the "norm" 

||A P ||(f) := ,/W Vfaxfdx, (8) 
y j Jn>o 

where the summation includes the variables in the tuple (<?,/), but does not include the 
variables A. 



D. The exact solutions used in the numerical experiments 

Only a few exact solutions of Einstein's vacuum equations possess both, high symme- 
tries and time dependence, or even better gravitational radiation. The asymptotically A3 
solutions do. They can be written as 

g = ^ e M (-dt 2 + dx 2 ) + l(t 2 + x 2 )(e w dy 2 + e~ w dz 2 ) (9a) 
\ft 2 + x 2 2 

n = ±(t 2 -x 2 ), (9b) 

where the functions M(t, x) and W(t, x) are solutions of certain differential equations (cf. 0). 
We restricted our analysis to two cases, M = and W = 0, the A3 solution on the one 
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hand, and M = — ^(t 2 + x 2 ) 2 and W = |(t 2 — x 2 ) on the other hand. The second case, 
unlike the first one, contains gravitational radiation. Both solutions behaved very similarly 
in our numerical experiments. The amount of gravitational wave content does not seem to 
be significant for the drift away from the constraint submanifold. For shortness we only 
present calculations done with the A3 solution. 

The solutions given in (|9|) are already extended beyond the two null-infinities J\ and J 2 
at t — —x and t = +x, the shaded region in FIG. [1] shows the physical part. In this repre- 
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sentation the point (0,0) represents future time-like infinity i + . The metric and curvature 
quantities diverge there. When approaching this point in our numerical evolution, the ab- 
solute value of several quantities must increase strongly, which implies increasing absolute 
errors. However, we confirmed that, even when approaching i + , our numerical scheme is 
second-order convergent, the difference between the numerical and the exact solution as well 
as the violation of constraints converges with a convergence rate of two. In the numerical 
experiments we start our calculation at t — — 1 and stop at t = —1/2. 
In other coordinates FIG. p] looks different, e.g. one can choose the coordinates such that 
the jTs are at constant x' values and that i + lies at a conformal time t' = oo. Since our 
findings are not influenced by this change of coordinates we refrain from a presentation. 



IV. THE QUALITY OF THE NUMERICAL SOLUTION: PARAMETER STUDY 

In this section we discuss the effect of the various parameters in the A-system @. In the 
case of the ID conformal field equations, the matrices B, C, D and E take values in R 14,14 , 
R 7 ' 14 , R 14,14 , and R 14 ' 14 . Since the parameter space is infinite, we had to restrict ourselves 
to exemplary cases, an exhaustive study by numerical means is impossible. 
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A. Overview 



Before going into detail, we discuss the main result of our numerical experiments. Al- 
though we were not able to find a suitable choice of parameters such that the constraint 
submanifold became an attractor for all times, we were able to improve the violation of con- 
straints up to a factor 5 compared to a numerical evolution of the unextended system. There 
are strong numerical indications that this improvement is obtained at the cost of a solution 
which is worse with respect to the pseudo difference norm. As an example, the dashed lines 
in FIG. |2] and |3] show the norms for the violation of constraints (0) and the pseudodiffer- 
ence (|8|) for a choice of parameters B = 3 • 1, C = 0,D = 0, and E = 1 in the A-system (by 
1 we denote the 14 x 14 unit matrix). The solid lines show the corresponding results for an 
evolution of the unextended system ([!]). FIG. ^| shows that the A-system initially attracts 
towards the constraint submanifold (until t = —0.85), before the violation increases with 
time. FIG. |3] shows that the pseudodifference norm at the same time increases more rapidly 
right from the beginning than in the numerical solution of the unextended system. 




FIG. 2: Constraint norm ||Ac||(i) for the A- 
system of run 7 (dashed line) in comparison 
to the unextended system (solid line). 



FIG. 3: Pseudodifference norm ||A-p||(i) for 
the A-system of run 7 (dashed line) in com- 
parison to the unextended system (solid line) . 



In TABLE 1, we present a summary of the performed numerical experiments. The num- 
bers in the parameter columns denote the value of the diagonal elements of the corresponding 
diagonal matrix. In the cases in which we studied a range of diagonal elements, we condense 
by combining to parameter intervals with similar behaviour. The observed development 
of the constraint norm ||Ac|| and the pseudodifference norm ||A-p|| is described using the 
following notation: 

— / +: The norm is smaller/greater than that of the unextended system in the whole domain 
of time integration. 

t / |: The norm is smaller/greater than that of the unextended system after small integra- 
tion times and greater /smaller at the end of the integration. 

The parameter choices for the run 1 were based on our expectations explained in section II. 
With this choice we were able to reduce the growth of the violation of constraints, but we 
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a) C = (Ci,C 2 ) e R 7,14 , where Ci,C 2 € R 7 ' 7 . 
/?) Ia 2 denotes a diagonal matrix with 1 + Af at 
7) all other elements vanish. 

5) x denotes a diagonal matrix with space depending diagonal elements. 



TABLE I: Results of the numerical experiments for various A-systems. A num- 
ber in the column parameter denotes a diagonal matrix with the number on 
the diagonal. For the notation used in the column norm, please refer to the 
text. 

neither made the constraint submanifold an attractor nor did we improve the pseudodiffer- 
ence at the end of the evolution. Since we observed that the variables A grew faster and 
to larger values than we had expected, we added terms proportional to A 2 (run 2-4), to 
increase the damping for non-vanishing As. We did not observe any significant change in 
the behaviour. This may have two reasons: Either the additional damping is too weak, a 
change in the parameters D and E alone is not sufficient, or these parameters are not the 
appropriate slots. To obtain a better understanding of the effect of the single parameters 
we performed the numerical experiments 5-19 whose results we are going to discuss in the 
following subsections. 
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B. Influence of the parameter B 

In the experiments 5-7, we studied the influence of the parameter B setting B propor- 
tional to the unit matrix, B = b 1, and varying b. The matrix E equals the unit matrix in 
all cases. The choice of the parameter b changes some of the characteristics of the A-system. 
Absolute values of b greater than 1 imply that the A-system has characteristic speeds larger 
than the speed of light. To avoid any influence of the grid boundary treatment on the 
numerical solution in the physical part of the grid — in the case |6| > 1 the outer grid 
boundary is no longer causally disconnected from the physical part of the grid — we moved 
the grid boundaries further out, something one can easily afford to do in a ID calculation 
with moderate values b. For our purposes this was sufficient and we, therefore, did not dis- 
cretise the analytic treatment of the initial boundary value problem for the A-system which 
would guarantee that no constraint violation is fed in from the grid boundary. For |6| > 1 
the maximally allowed time step is smaller than that of the unextended system due to the 
Courant-Friedrichs-Lewy condition. To compare numerical experiments for systems with 
different values b, we use the same, most restrictive time step for all runs. 

For small values of b we do not observe a significant change in the behaviour of the A- 
system compared to the unextended system, as can be seen in the FIG. f| and §. In these 




FIG. 4: Constraint norm ||Ac||(i) for the runs 
5 with b = — 1 (dashed line) and 6 (dotted 
line) in comparison to the unextended system 
(solid line). 



FIG. 5: Pseudodifference norm ||A-p||(i) for 
the runs 5 with b = — 1 (dashed line) and 6 
(dotted line) in comparison to the unextended 
system (solid line). 



figures we plot the constraint norm ([/]) and the pseudodifference norm (§) for b = — 1 (dashed 
line) and for b = 1.2 (dotted line) in addition to the unextended system (solid line). The 
curves are very similar, as are all other curves from the runs summarised as run number 5 
in TABLE |. There is a slightly reduced growth in the constraint norm and an increase in 
the pseudodifference norm at the end of the evolution. 

For a value of b significantly greater than 1 (run 7 and FIG. |2] and |3]) the constraint 
norm decreases initially and is always smaller than in the unextended system. In contrast, 
as already mentioned, the pseudodifference norm is always larger than in the unextended 
system. 
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C. Influence of the parameter E 

In runs 8-14 we studied the role of the parameter E, fixing the value of B to B = 1. A 
non-zero value of B is indeed necessary if A is to measure the violation of the constraints. 
We choose the value B = 1, as the characteristic speeds determined by B then agree with 
speed of light which seems to be a natural choice for the field equations of general relativity. 
In addition, the results in the previous subsection were rather insensitive to the exact value 
of the parameter B in this part of the parameter space. Changing the value of diagonal 
elements for a constant diagonal matrix E in experiments 8-11, we found the best results 
for the violation of the constraints and the pseudo difference for a value of E = —3 1 (FIG. |6| 
and FIG. |7]), where the variables A are amplified by the EA-term in the A-system. With this 
choice of parameters, the pseudodifference norm can be slightly improved during the whole 
numerical integration up to t = —0.5, but we stress that this improvement is not significant. 
We checked the results for this run after a longer integration time (up to t = —0.3). We then 
found a worse pseudodifference norm compared to the unextended system. In runs 12-14 we 
put only single diagonal elements of E to -3, (E)n and (E) 88 , which influence directly those 
two constraints which are most vehemently violated. These constraints are the constraints 
for the quantities hu and En. Only affecting both constraints with our choice of parameters 
can improve the numerical evolution of the constraints and the pseudodifference compared 
to the unextended system. However, in comparison to run 10, the results are worse. 
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FIG. 6: Constraint norm ||Ac||(t) for run 
10 in comparison to the unextended system 
(solid line). 
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FIG. 7: Pseudodifference norm ||A-p||(£) for 
run 10 in comparison to the unextended sys- 
tem (solid line). 



D. Influence of the parameter D 

With the experiment 15, we studied the influence of a constant, diagonal matrix D 
together with a matrix B = 1. In all runs with non-vanishing D the violation of constraints 
could be improved (see e.g. FIG. § for d — 10) at the cost of a worse pseudodifference. 
The constraint and the pseudodifference norm are identical for runs with D = d 1 and 
D = —d I. This property results from a symmetry of all time evolution variables in (|3]) 
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under the simultaneous transition D — > — D, x — > — x for this specific choice of parameters. 
For the A3-solution, all evolution variables at a fixed time are even/odd functions on the 
space coordinate x. For a diagonal, constant matrix D, the term DA couples to an even/odd 
function /, the corresponding Aj, which has, according to ([k]), the opposite symmetry, as it 
measures the corresponding constraint, which involves one space derivative of f\. 
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FIG. 8: Constraint norm ||Ac||(i) for run 15 
with d = 10 in comparison to the unextended 
system (solid line). 



FIG. 9: Pseudo-difference norm ||A-p||(i) for 
run 19 in comparison to the unextended sys- 
tem (solid line). 



E. Influence of the parameter C 

In the experiments with a non-vanishing matrix C (run 16) the violation of the constraints 
was only improved at the cost of a worse pseudodifference. As the result for the constraint 
norm is similar to that in FIG. [8], and the result for the pseudodifference norm is qualitatively 
given by FIG. [3], we refrain from presenting figures for run 16. 

F. Influence of the parameters D and E 

Motivated by the good results for the pseudodifference in run 10 and for the constraints 
in run 15, we studied the correspondence of non- vanishing, diagonal parameters D and E for 
a parameter B = 1 in runs 17-19. Using the same parameters B, C and E as in run 10, but 
now with a non-vanishing D, we could not improve the numerical solution in run 17. Hence, 
in runs 18-19, we again stuck to our original choice of E = 1. In run 18, the pseudodifference 
first increased, before approaching the curve of the unextended system at integration times 
of about t = —0.5 (FIG. ||). In experiment 19, we chose the matrix D = 5x 1 in order to 
keep the symmetries of the A3 solution. Again, the constraints were only improved at the 
cost of a worse solution. 
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V. CONCLUSION AND OUTLOOK 

The freedom in extending a system of evolution equations with constraints to a A-system 
is huge. For the conformal field equations of general relativity, we have explored the effect of 
what we thought are natural choices in the A-system, analysing the influence on the quality 
of the new system's numerical solution. We found that we were able to significantly reduce 
the violation of the constraints, but we also found that this improvement did not imply a 
smaller numerical error. 

Furthermore, we found that the significant improvement in the violation of the constraint 
did not prevent the solution from eventually running away from the constraint submanifold, 
i.e. the A-systems used do not inherit the property of asymptotic stability from the linear 
system. Similar experience with the semi-linear SU (2)- Yang-Mills equations [IT| and the 



study of a simplified model system suggest that a more balanced choice of parameters 



of the A-system is needed to achieve an attractive force towards the constraint submanifold 
for all times. Recent analytic results by Heinz-Otto Kreiss and Peter Hiibner give sufficient 
conditions for asymptotic stability. As we found a strong correlation between a smaller 
violation of the constraints and a worse numerical solution in the A-system we suspect, 
however, that asymptotic stability does not necessarily imply a smaller numerical error. 
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